Skip to content

Add CST parity functions and optimize FRD workflows#22

Merged
jamestjsp merged 8 commits intomainfrom
feat/matlab-cst-parity
Apr 3, 2026
Merged

Add CST parity functions and optimize FRD workflows#22
jamestjsp merged 8 commits intomainfrom
feat/matlab-cst-parity

Conversation

@jamestjsp
Copy link
Copy Markdown
Owner

Summary

  • add the MATLAB Control System Toolbox parity surface and the accompanying benchmark suite
  • fix follow-up correctness issues in Lsim, FRDMargin, and modal eigenvalue classification
  • optimize System.FRD, FRD.Sigma, FRDFeedback, and related FRD benchmark paths, and rename bench_new_test.go to benchmark_hotspots_test.go

Testing

  • go test -count=1 ./...
  • go test -count=1 -race ./...
  • go test ./... -run 'TestCanon'\n- go test -run '^$' -bench 'BenchmarkSystemFRD_MIMO_2000|BenchmarkCanon_Modal_N50' -count=1 ./...\n

jamestjsp and others added 8 commits April 2, 2026 20:51
Phase 1 — Quick wins:
  Norm, Covar, Lsim, Augstate, SS2SS, Xperm, Inv, Pzmap, Isproper

Phase 2 — PID controllers:
  NewPID (parallel form), NewPIDStd (standard form), PID2 (2-DOF),
  Pidtune (auto-tuning: P/I/PI/PD/PID/PIDF via open-loop shaping)

Phase 3 — Analysis & canonical forms:
  Canon (modal + companion), Loopsens, Rss/Drss, Ssbal, Prescale

Phase 4 — Frequency response data model:
  FRD type, FRDSeries/FRDParallel/FRDFeedback, FRD Nyquist/Sigma/Margin

Phase 5 — Structural decomposition:
  Sminreal, Stabsep, Modsep (shared eigdecomp core with Schur fallback)

Also fixes README inaccuracies (Norm, Lsim, InitialCondition).

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
1. FRD.Sigma() MIMO: real block form doubles each SV; take every-other
   value (vals[2*s]) instead of the first min(p,m). diag(2,1) now
   correctly returns [2,1] not [2,2].

2. Delay metadata: Inv, Canon, Ssbal, Prescale, Sminreal, Stabsep, and
   Modsep now reject systems with delays (HasDelay) instead of silently
   dropping Delay/InputDelay/OutputDelay/LFT metadata.

3. Loopsens(P, C): changed from Loopsens(L) to accept plant and
   controller separately. Computes true input-side sensitivities
   Si=(I+CP)^{-1}, Ti=CP(I+CP)^{-1} via C*P loop, which differ from
   output-side So/To for non-commutative MIMO. Added regression test
   verifying Si != So.

4. FRDMargin: rewrote to use Bode()'s unwrapped phase instead of raw
   cmplx.Phase, fixing missed crossings at ±180° wrap. Fixed WgFreq/
   WpFreq field assignment to match existing Margin() convention.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
FRDMargin: collect all crossings then apply same two-pass selection as
Margin() — prefer smallest positive margin, fall back to any finite
value. K=10/(s+1)^3 now returns GM≈-1.94dB matching Margin().

Lsim: validate that t vector is uniformly spaced (tol 1e-6 relative)
before discretizing with dt=t[1]-t[0]; reject non-uniform grids with
a clear error instead of silently producing wrong output.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Adds check in the discrete branch: abs(sys.Dt - dt)/sys.Dt <= 1e-6,
so a discrete sys with Dt=0.1 fed a t grid at 0.05 spacing now errors
instead of silently simulating with the wrong sample period.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
FRD stack (System.FRD, Sigma, Margin, Series, Parallel, Feedback):
  SISO/MIMO nw=200/2000/10000, large n=50 MIMO, 5x5 MIMO,
  negative-margin FRDMargin case

Time response (Step, Impulse, Lsim):
  SISO n=2/10/50, MIMO n=10, Lsim 1e4 steps, discrete Lsim

Loopsens: SISO n=2/10, MIMO n=10

Transform/decomposition (Canon, Stabsep, Modsep, Ssbal, Prescale,
Sminreal, Inv, Covar): n=2/10/50/100 sweeps, SISO + MIMO

Pidtune: PI, PID, PIDF on 4-state plant

All fixtures use non-symmetric A (sub+superdiagonal coupling).

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
System.FRD: single-alloc contiguous response storage via
newFRDResponseStorage; bulk copy from FreqResponse.Data.

FRD.Sigma: pre-allocate SVD workspace (block matrix, work array,
singular values) once via newFRDSVDWorkspace, reuse across all
frequency points; call impl.Dgesvd directly instead of mat.SVD.

FRDSeries/Parallel: write into pre-allocated flat data slice via
cMulNestedInto/cAddNestedInto instead of per-frequency [][]complex128.

FRDFeedback: pre-allocate all scratch buffers (G, K, KG, I+KG, inv,
aug) in newFRDFeedbackWorkspace; flat complex128 multiply/invert via
cMulInto/cInvertInto; zero per-iteration allocations.

benchstat (geomean): -14% time, -22% memory
  SystemFRD SISO nw=10k:     -50% time, -99.7% allocs
  FRDSigma MIMO nw=10k:      -14% time, -93% memory
  FRDFeedback SISO nw=2k:    -74% time, -58% memory
  FRDFeedback MIMO nw=10k:   -54% time, -81% memory
  FRDSeries MIMO nw=10k:     -39% time, -99.9% allocs

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Absolute 1e-10 tolerance breaks for eigenvalues with magnitude >>1
where machine eps on the imaginary part exceeds the threshold. Use
1e-10*max(1, |lambda|) consistent with eigdecomp.go's isConjugate.

Co-Authored-By: Claude Opus 4.6 (1M context) <noreply@anthropic.com>
Search for a conjugate partner before collapsing a small-imaginary eigenvalue to a real modal block, add a regression test for the large-real/small-imag pair case, and rename bench_new_test.go to benchmark_hotspots_test.go.

Co-authored-by: Copilot <223556219+Copilot@users.noreply.github.com>
@jamestjsp jamestjsp merged commit 0800e07 into main Apr 3, 2026
2 checks passed
@jamestjsp jamestjsp deleted the feat/matlab-cst-parity branch April 3, 2026 14:55
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant